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ABSTRACT 

We have used maximum-likelihood estimation to determine the 
quadrupole amplitude Q r ms-ps and the spectral index n of the density 
fluctuation power spectrum at recombination from the COBE DMR data. 
We find a strong correlation between the two parameters of the form 
Qrms-ps = (15.7 ± 2.6) exp[0. 46(1 — n)) /xK for fixed n. Our result is slightly 
smaller than and has a smaller statistical uncertainty than the 1992 estimate of 
Smoot et al. 



Subject headings: cosmic microwave background — cosmology: observations 
large-scale stucture of universe 



1. Introduction 



The search for the cosmic microwave background radiation (CMBR) anisotropies 
finally yielded positive results with the detection made with the Differential Microwave 
Radiometers (DMRs) on the Cosmic Background Explorer (COBE) satellite (Smoot et al. 
1992; Bennett et al. 1992; Wright et al. 1992). These authors analyzed sky maps based on 
the first year of COBE data and presented quantitative results for the quadrupole moment, 
the angular correlation function, and power spectrum parameters characterizing the large 
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angular scale fluctuations. In simple cosmological models the CMBR anisotropy is related 
directly to the gravitational potential fluctuations present during recombination on large 
length scales ( [Sachs fc Wolle 1967| ). The COBE data thus enable quantitative testing of 
cosmological theories, free from the complications of galaxy formation. 

Perhaps the most important measure of CMBR anisotropy is provided by the angular 
correlation function of the temperature anisotropy AT, C(a) = (AT(ra 1 )AT(ri2))- The 
average is taken over the sky for all pairs of directions separated by angle a = cos^ 1 (?i! • n^). 
Smoot et al. (1992, Fig. 3) and Wright et al. (1992, Fig. 2) presented estimates for C(a) 
with the monopole and dipole contributions removed; Smoot et al. (but not Wright et al.) 
also removed the quadrupole. 

Smoot et al. (1992) used a least-squares fit to the measured values C(a>i) to estimate 
the amplitude (represented by the rms CMBR quadrupole, Qrms-ps) and the spectral 
index n of the density power spectrum, obtaining n = 1.1 ± 0.5 and Q rm s-ps = 17 ± 5 /iK 
for n — 1, corresponding to the scale-invariant Peebles-Harrison-Zel'dovich spectrum. In 
order to make these estimates one needs the covariance matrix for C(a.j). There are two 
sources of uncertainty contributing to this covariance matrix: measurement errors (chiefly 
receiver noise) and cosmic variance (the intrinsic statistical fluctuations expected because 
the CMBR temperature is measured on a surface of finite extent). The cosmic variance, and 
therefore the full covariance matrix of the measurements, depends on the power spectrum 
parameters one is trying to determine. Smoot et al. presented least-squares estimates with 
and without the cosmic variance. The estimated values changed very little but the x 2 
decreased from 79 to 53 for 68 degrees of freedom when the cosmic variance was included. 
We find that the inclusion of cosmic variance is even more important for the correlation 
function of Wright et al. (1992, Fig. 2), as it causes \ 2 to decrease from 155 to 52 for 64 
degrees of freedom if n — 1 and Q rms „ps = 17 /iK. 

Recently, Scaramella & Vittorio (1993) emphasized the importance of cosmic variance. 
Including this, they reanalyzed the angular correlation function of Smoot et al. (1992) using 
X 2 minimization and Monte Carlo simulations. They concluded that the best-fit amplitude 
is Qrms-ps — (14.5 ± 1.7) (1 ± 0.06) /iK for n = 1, where the first error bar is due to cosmic 
variance and the second one to measurement uncertainties. However, they did not simulate 
the actual sky sampling and data reduction procedure applied to the data. 

There are several shortcomings of the previous statistical analyses. First, the least- 
squares method is inappropriate when the covariance matrix depends on the parameters one 
is trying to estimate. In this case least-squares estimators are often biased and in general 
they do not have minimum variance. The standard least squares error estimates are also 
biased, and the sum of squares of residuals does not have a x 2 distribution. Correlations 
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between different lag angles a, are not taken into account because the least-squares method 
(and also the minimum \ 2 method) uses only the diagonal part of the covariance matrix 
of C(ttj). As we show below, this covariance matrix has three contributions, two of which 
involve cosmic variance and are nondiagonal. Although the biases can be corrected using 
Monte Carlo simulations, least-squares estimators do not make the most efficient use of the 
data because they are not minimum-variance. 

Wright et al. (1993) reexamined the COBE data using the rms variance on the 10° 
scale and the Boughn-Cottingham statistic (Boughn et al. 1992). They confirm their earlier 
results for the quadrupole amplitude, concluding that the amplitude is consistent with 
Qrms-ps = 17 ± 5 /iK for n = 1. 

We have chosen instead to estimate the power spectrum parameters using the 
maximum-likelihood method. This method has the advantage of providing estimates that 
are asymptotically unbiased and of minimum variance (Eadie 1971). This Letter presents 
our approximate maximum-likelihood determination of Qrms-ps an d n using realistic Monte 
Carlo simulations of the COBE sampling and data reduction procedures. We will show 
that the uncertainties of these two parameters are strongly correlated and that the optimal 
value of Qrms-ps(^) for the correlation functions of Smoot et al. (1992) and Wright et al. 
(1992) is slightly lower and has smaller uncertainty than estimated by those authors. 



2. Statistical Method 



The measured correlation function C(ai) was given in m (= 70 for Smoot et al. 1992 
and 65 for Wright et al. 1992) equidistant points between 0° and 180°. Approximating 
the m-dimensional distribution as a multivariate normal distribution one can write the 
likelihood function as 

exp [-|(AC) T M- 1 (AC)] 

L(Q rms -PS, n) = [(27r)mdet(M)] i/ 2 • C 1 ) 

Here (AC) T and (AC) are m-dimensional row and column vectors with entries 
(AC), = C(ai) - (C(ai)), while M = ((AC) (AC) T ) is an m x m matrix. The angle 
brackets here denote averages over both measurement errors and the statistical ensemble of 
temperature fluctuation fields for given Qrms-ps and n. Both (C(a)) and M depend on these 
parameters. Maximum-likelihood estimates Q rm s-ps and n are obtained by maximizing 
L(Qrms-PS) n ) 5 keeping C{a.j) fixed at the measured values. The curvature of the likelihood 
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function around the maximum provides an estimate of the covariance matrix of errors of 
the estimated values (Eadie 1971). 

In assuming a multivariate normal distribution we are maximizing the wrong likelihood 
function because the actual distribution of C(oti) involves products of normally distributed 
variables (the raw temperature measurements) and in general is not normally distributed. 
However, the assumption of a normal distribution is not crucial to the success of our 
method when there are small departures from normality. Although the asymptotic property 
of the smallest variance among all estimation methods is no longer valid, the increase in 
variance should be small. One could appeal to the central limit theorem and assume that 
the average over many products of pairs is nearly normally distributed even when the 
individual products themselves are not. We have tested this possibility by making Monte 
Carlo realizations of pixel maps and calculating from them the distributions of C{a.j). 
They are close to normal away from the tails. The extended tails increase the variance of 
the estimated values, but the effect is small and does not significantly affect the overall 
efficiency of the maximum-likelihood method. We will use Monte Carlo simulations with 
the correct distributions to estimate the bias and variance of Q rms -ps and n. 

The DMR instrument measured, for two independent channels at each of three 
frequencies, the differences in CMBR temperature Tj — Tj between two directions Hi and n,j 
separated by 60° (Smoot et al. 1990). These differences were then fitted to give estimates 
of ATj (with the monopole, dipole, and possibly quadrupole contributions removed) at 
6144 points of an equal area sky map (Torres et al. 1989; Janssen & Gulkis 1992). The 
correlation function was then calculated by Smoot et al. (1992) and Wright et al. (1992) 
using 

, ) = J2i N t Ej Nj m ATj 5(n t ■ g - cos a k ) 

{ak> Ei N i Ej Nj 5(ni ■ nj - cos a k ) ' U 

where TVj is the number of measurements in each pixel i and the 5-function indicates that the 
sum is performed over all different pairs ATj and ATj such that cos _1 (rTj -nj) is in a given bin 
k of a. The weighting by iVj gives a minimum- variance estimate for C{a,k). The temperature 
anisotropies ATj and ATj can be drawn either from the same map (autocorrelation) or 
from two different maps (cross-correlation). Smoot et al. cross-correlated the 53A+B and 
90A+B GHz channels, while Wright et al. combined cross-correlations of 53A, 53B, and 
90B channels in a manner equivalent to the autocorrelation of a single weighted map. 

Evaluating the mean and covariance matrix of C(ak) requires knowing the TVj and 
the covariance matrix of temperature measurement errors. These quantities depend in a 
complex way on the COBE sky scan pattern and on the detailed properties of the DMR 
instrumentation (Boggess et al. 1992; Smoot et al. 1990). We used a simulation program, 
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kindly provided by Ned Wright, to simulate the spacecraft operation for the first year of 
operation, including the correct orbit, spacecraft spin, tracking of the Sun, Moon, and 
planets, with data rejection if the instrument pointed too close to the Earth and Moon, 
etc. (Smoot et al. 1992). The simulated measurements were gathered into 6144 equal area 
pixels using the quad-cube routines provide to us by Wright. 

In addition, we calculated the covariance matrix of temperature measurement errors 
for the sky maps. The main source of measurement error is receiver noise (Smoot et al. 
1990). Because the raw DMR measurements are temperature differences for two beams 
separated by 60°, the errors in the temperature maps (obtained by fitting to the differences) 
are correlated. We estimated that the off-diagonal elements of the covariance matrix are 
nearly all much smaller than 1% of the diagonal elements (rising to 6% for a few elements) 
so that to a good approximation one can safely neglect the noise correlations for different 
pixels. One can then write ATj = AT! + where ATj is the measured value for pixel i, 
AT the true value and q the noise contribution. The pixel noise is normally distributed 
with (ej) = and (titj) = (c 2 /Ni)5ij, where a is the noise contribution from a single 
measurement (Janssen & Gulkis 1992). 

The true values AT are also stochastic variables in theories of large-scale structure. 
The ensemble averages are (A?; ) = and (ATf AT ) = C°(a k ), where C°{a k ) is 
a theoretical correlation function (including beam smearing and pixelization). This 
function is most conveniently characterized by its expansion in Legendre polynomials, 
C°(a k ) = J2i Ci Gf Pi(cosctk), where Gi is the window function of the DMR beam, for 
which we used the values given by Wright et al. (1993) with a slight correction for beam 
smearing and pixelization. The angular power spectrum on large angular scales for primeval 
adiabatic density fluctuations with = 1 is (Bond & Efstathiou 1987) 

(2/ + l)r[/ + (n-l)/2]r[(9-n)/2] 
' 5 r[/ + (5-n)/2]T[(3 + n)/2] ^ rms - ps ' { ) 

This expression is accurate for the angular scales probed by COBE. To test it we replaced 
it with the more accurate angular correlation function obtained with a full integration of 
the coupled Boltzmann and Einstein equations for n — 1 by Bond & Efstathiou (1987) and 
found negligible change in our estimates. 

Averaging over measurement errors and an ensemble of true sky maps we can now 
calculate the mean and covariance matrix of C(a k ) neglecting the fitting and removal 
of low-order multipoles. For the cross-correlation case we get (C(ctk)) = C°(a k ); for 
autocorrelations there is an additional term at a k = due to noise. The covariance matrix 
for the cross-correlation case is the sum of three terms, M = NN + SN + SS, which we denote 
as noise-noise (NN), signal-noise (SN) and signal-signal (SS) terms. The NN term depends 
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only on the measurement errors, the SN term scales as Q^ms-ps an d the measurement 
variance, and SS scales as Qtms-ps- The contributions to the matrix elements M klk2 are 
given by the following expressions: 



NN = ?A^hh £ Nn £ Nn 5{Hii . n n _ cos aki ) , (4 ) 
iVtot h h 



SN = -!— N h Y N h 5{n h ■ n h - cos a kl ) 

iVtot h jl 



X 



°i Z ^72 <K™»i ' ™?2 - cos «fe 2 ) C°(cos 1 [ra^ • n i2 ]) 



J'2 



X 



+°b S ^2 <K™* 2 ' «ji ~ cos a fc 2 ) C°(cos 1 [rii 1 • n i2 ]) I , (5) 
SS = — !— Z ^"n Z % S ^2 Z % ^ (wii • ™n - cos a fcl ) <5 (n i2 • n i2 - cos a fc2 ) 

tot jj jj j 2 j 2 

{c°(cos _1 K • ^ 2 ]) C7°(cos- 1 [^- 1 • n n \) + C°(cos _1 K • n j2 ] ) C°(cos _1 [^i • ^ 2 ])} , (6) 
where 

A^tot = Y N h 2 % 2 ^2 Z ^J'a 5 (™n ' «ji ~ cos ) 5 (^2 ' "ja ~ cos «fc 2 ) 
U jl * 2 32 

and the indices range over all map pixels. Pixels labeled with i correspond to map A, 
for which the measurement variance of ATj is <J A /Ni, while j corresponds to map B. For 
cross-correlations a a ^ &b in general. For autocorrelations one sets <ja = &B and the NN 
and SN terms are increased by a factor of 2. 

We see that the NN term involves double summation over sky maps, the SN term 
involves triple summation, and the SS term quadruple summation. Even after the galactic 
lattitude cuts made by Smoot et al. (1992, \b\ > 20°) and Wright et al. (1992, \b\ > 30°), 
one is still left with several thousand pixels. While the NN and SN terms can be summed 
exactly, the direct calculation of SS becomes computationally too expensive. Instead, we 
evaluated it using Monte Carlo simulations. We generated 10,000 maps for a gaussian 
random field AT(n) on the sphere having each theoretical angular power spectrum (i.e., 
value of n) to be tested. For each realization the angular correlation function was measured 
using equation (2) and the ensemble average was made over the 10,000 samples. Note that 
by adding the noise to the Monte Carlo samples one could similarily calculate the total 
covariance matrix (as we did for testing). The advantage of dividing the whole covariance 
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matrix into three terms is that once we calculate the expression for one value of Q^ms-ps we 
can simply scale it to obtain the results for all different values of <3? ms _pg for a given n. The 
n-dependence of C°(a) is sufficiently smooth so that we interpolated the matrix elements 
of SN and SS evaluated on a grid of values of n. To test the whole procedure including 
the values of Ni we have compared the NN term with the measurement errors for C(ai) 
obtained by Smoot et al. (1992, Fig. 3) and Wright et al. (1992, Fig. 2). For both data 
sets our results agree with the correct values within a few percent. 

Given the full covariance matrix it is now straightforward to obtain maximum-likelihood 
estimates Qrms-ps and n for a given set of data C(otk)- The covariance matrix of errors in 
the parameters may be estimated in the usual way by taking AlnL = 0.5 for one standard 
deviation. However, we should not trust these asymptotic results because of the small 
numbers of independent data points given the COBE beam and the intrinsic correlations as 
well as our assumption of a normal likelihood function. 

There is another reason why our estimator may give biased results. Our procedure so 
far does not correctly simulate the data reduction procedure used by Smoot et al. (1992) 
and Wright et al. (1992) because we have not accounted for the fitting and removal of 
low-order multipole moments from the maps before the angular correlation function is 
computed. If the sky sampling were uniform and complete, this could be accounted for 
simply by limiting the range of / used to compute C°(ak)- However, incomplete sky coverage 
couples different multipoles so that the actual (C(ak)) and covariance matrix differ from 
what we give above, resulting in a bias in our maximum-likelihood estimator. 

To correct for this and other biases and to determine the variance of our estimator we 
resort again to Monte Carlo simulations. We generate 5000 random sky maps including 
signal and noise using the estimated values of Q rms -ps and n for the COBE data as the 
input parameters. We fit and remove monopole, dipole, and (optionally) quadrupole using 
the correct galactic latitude cut and then compute the angular correlation function. These 
samples are used as the input data in the maximum-likelihood estimation described above. 
The distribution of Monte Carlo estimates around the true value yields the variance and 
bias of the estimate. We use this bias estimate to correct our results given below. 



3. Results 

Using the Smoot et al. (1992, Fig. 3) data we obtain maximum-likelihood estimates, 
before bias correction, of n = 1.2 and Qrms-ps = 12.2 /iK. Using the Wright et al. (1992, 



- 8- 



Fig. 2) data the corresponding values are n = 0.9 and Qrms-ps = 13.9 /iK. As indicated 
above, our estimator is expected to be biased. To compensate for this bias we analyze 
Monte Carlo simulations of the data that are made with different choices for the true 
(n, Qrms-ps) and we try to determine those parameters for which the mean estimated values 
equal the ones we obtain from the real data. We find that the bias in n is less than 0.1, but 
the bias in Q rms -ps is significant. In the Smoot et al. case our bias for true values n = 1 and 
Qrms-ps — 15.0 /iK is AQrms-PS = —2.6 fiK (the mean estimate is (Qrms-ps) — 12.4 /iK) 
while in the Wright et al. case it is AQrms-PS = —2.2. The bias is larger for the Smoot 
et al. analysis because of the additional quadrupole subtraction applied to the data. The 
bias is only weakly dependent on n and Qrms-ps- Thus, the bias-corrected estimates are 
(™,<2rms-ps) = (1-2, 14.8 /iK) for Smoot et al. and (0.9, 16.1 /iK) for Wright et al. 

In general one is interested in the amplitude of fluctuations for a fixed value of n. 
Assuming the scale-invariant slope n = 1 and combining the two data sets we obtain 
Qrms-ps = 15.7 ± 2.6 ^K, where the uncertainty is taken from our fit to the Wright et al. 
sample. There is a strong anticorrelation between our estimates of n and Q rms -ps and they 
cannot be independently determined with high precision. We find that the approximate 
relation between the two parameters is of the form 

Q rm s-ps = (15.7 ± 2.6) exp[0.46(l - n)\ /iK . (7) 

Our mean value is slightly higher than that obtained by Adams et al. (1993) using the 
<7 s k y (10°) normalization, Q rm s-ps = 15(1 ± 0.2) exp[0.31(l — n\) /iK. It is also higher than 
that obtained by Scaramella & Vittorio (1993). The main reason for these differences is that 
the fitting and removal of the dipole and quadrupole applied to the real data also subtracted 
some of the higher-order multipole moments because of nonuniform sky coverage. 

Our results for n — 1.0 agree within the errors with the results reported by Smoot et 
al. (Qrms-ps = 16.7 ±4.7 //K), with a slightly lower amplitude and a smaller error bar. The 
change in the amplitude, when combined with the smaller error bar, may not be enough 
to significantly improve the consistency of the COBE results with the low upper limit on 
smaller angular scales obtained recently by Gaier et al. (1992). 

This work was supported by grants NSF AST90-01762 and NASA NAGW-2807. We 
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form, and we acknowledge helpful discussions with C. Lineweaver and C. Bennett. We are 
especially thankful for the comments and assistance given by the referee, Ned Wright. 
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